no_sp_model <- organ_sp %>% 
  lmer(eligible_population_enrolled ~ total_days + opo + (1 | county), data = ., 
       REML = FALSE)

no_sp_model %>% 
  broom::tidy() %>% 
  filter(term == "total_days") %>% 
  mutate(estimate_year = estimate * 365, 
         se_year = std.error * 365, 
         ci_low = estimate_year - 1.96*se_year, 
         ci_upper = estimate_year + 1.96*se_year) %>% 
  mutate(term = ifelse(term == "total_days", "Time (year)", term)) %>% 
  select(term, estimate_year, se_year, ci_low, ci_upper) %>% 
  rename("Term" = term, 
         "$\\beta$" = estimate_year, 
         "Standard error" = se_year, 
         "Lower bound" = ci_low, 
         "Upper bound" = ci_upper) %>% 
  knitr::kable(digits = 3, escape = F) %>% 
  kable_styling(full_width = F, bootstrap_options = c("hover", "striped"))

Term	β	Standard error	Lower bound	Upper bound
Time (year)	2.675	0.008	2.66	2.69

organ_sp %>% 
  ggplot(aes(x = total_days, y = eligible_population_enrolled)) + 
  geom_jitter()

add_predictions(organ_sp, no_sp_model) %>% 
  ggplot(aes(x = total_days, y = pred, color = opo, type = county)) + 
  geom_line() + 
  geom_jitter(aes(x = total_days, y = eligible_population_enrolled), color = "black", alpha = 0.1) +
  viridis::scale_color_viridis(discrete = TRUE) + 
  theme(legend.position = "bottom")

sp_12_model <- organ_sp %>% 
  lmer(eligible_population_enrolled ~ total_days + year_sp_2012 + (1 | county), data = .)

sp_12_model %>% 
  broom::tidy() %>% 
  filter(term %in% c("total_days", "year_sp_2012")) %>% 
  mutate(estimate_year = estimate * 365, 
         se_year = std.error * 365, 
         ci_low = estimate_year - 1.96*se_year, 
         ci_upper = estimate_year + 1.96*se_year) %>% 
  mutate(term = ifelse(term == "total_days", "Time (year)", term), 
         term = ifelse(term == "year_sp_2012", "2012 spline", term)) %>% 
  select(term, estimate_year, se_year, ci_low, ci_upper) %>% 
  rename("Term" = term, 
         "$\\beta$" = estimate_year, 
         "Standard error" = se_year, 
         "Lower bound" = ci_low, 
         "Upper bound" = ci_upper) %>% 
  knitr::kable(digits = 3) %>% 
  kable_styling(full_width = F, bootstrap_options = c("hover", "striped"))

Term	β	Standard error	Lower bound	Upper bound
Time (year)	2.883	0.022	2.840	2.926
2012 spline	-0.326	0.032	-0.389	-0.262

add_predictions(organ_sp, sp_12_model) %>% 
  ggplot(aes(x = total_days, y = pred, color = opo, type = county)) + 
  geom_line() +
  viridis::scale_color_viridis(discrete = TRUE) + 
  theme(legend.position = "bottom")

sp_16_model <- organ_sp %>% 
  lmer(eligible_population_enrolled ~ total_days + year_sp_2016 + opo + (1 | county), data = .) 

sp_16_model %>% 
  broom::tidy() %>% 
  filter(term %in% c("total_days", "year_sp_2016")) %>% 
  mutate(estimate_year = estimate * 365, 
         se_year = std.error * 365, 
         ci_low = estimate_year - 1.96*se_year, 
         ci_upper = estimate_year + 1.96*se_year) %>% 
  mutate(term = ifelse(term == "total_days", "Time (year)", term), 
         term = ifelse(term == "year_sp_2016", "2016 spline", term)) %>% 
  select(term, estimate_year, se_year, ci_low, ci_upper) %>% 
  rename("Term" = term, 
         "$\\beta$" = estimate_year, 
         "Standard error" = se_year, 
         "Lower bound" = ci_low, 
         "Upper bound" = ci_upper) %>% 
  knitr::kable(digits = 3) %>% 
  kable_styling(full_width = F, bootstrap_options = c("hover", "striped"))

Term	β	Standard error	Lower bound	Upper bound
Time (year)	2.563	0.01	2.542	2.583
2016 spline	0.789	0.05	0.692	0.886

add_predictions(organ_sp, sp_16_model) %>% 
  ggplot(aes(x = total_days, y = pred, color = opo, type = county)) + 
  geom_line() + 
  viridis::scale_color_viridis(discrete = TRUE) + 
  theme(legend.position = "bottom")

sp_17_model <- organ_sp %>% 
  lmer(eligible_population_enrolled ~ total_days + year_sp_2017 + opo + (1 | county), data = ., 
       REML = FALSE)

sp_17_model %>% 
  broom::tidy() %>% 
  filter(term %in% c("total_days", "year_sp_2017")) %>% 
  mutate(estimate_year = estimate * 365, 
         se_year = std.error * 365, 
         ci_low = estimate_year - 1.96*se_year, 
         ci_upper = estimate_year + 1.96*se_year) %>% 
  mutate(term = ifelse(term == "total_days", "Time (year)", term), 
         term = ifelse(term == "year_sp_2017", "2017 spline", term)) %>% 
  select(term, estimate_year, se_year, ci_low, ci_upper) %>% 
  rename("Term" = term, 
         "$\\beta$" = estimate_year, 
         "Standard error" = se_year, 
         "Lower bound" = ci_low, 
         "Upper bound" = ci_upper) %>% 
  knitr::kable(digits = 3) %>% 
  kable_styling(full_width = F, bootstrap_options = c("hover", "striped"))

Term	β	Standard error	Lower bound	Upper bound
Time (year)	2.561	0.009	2.543	2.579
2017 spline	1.700	0.078	1.548	1.852

add_predictions(organ_sp, sp_17_model) %>% 
  ggplot(aes(x = total_days, y = pred, color = opo, type = county)) + 
  geom_line() +
  viridis::scale_color_viridis(discrete = TRUE) + 
  theme(legend.position = "bottom")

all_sp_model <- organ_sp %>% 
  lmer(eligible_population_enrolled ~ total_days + year_sp_2012 + year_sp_2016 + year_sp_2017 + opo + 
         (1 | county), data = ., REML = FALSE)

all_sp_model %>% 
  broom::tidy() %>% 
  filter(term %in% c("total_days", "year_sp_2012", "year_sp_2016", "year_sp_2017")) %>% 
  mutate(estimate_year = estimate * 365, 
         se_year = std.error * 365, 
         ci_low = estimate_year - 1.96*se_year, 
         ci_upper = estimate_year + 1.96*se_year) %>% 
  mutate(term = ifelse(term == "total_days", "Time (year)", term), 
         term = ifelse(term == "year_sp_2012", "2012 spline", term), 
         term = ifelse(term == "year_sp_2016", "2016 spline", term), 
         term = ifelse(term == "year_sp_2017", "2017 spline", term)) %>% 
  select(term, estimate_year, se_year, ci_low, ci_upper) %>% 
  rename("Term" = term, 
         "$\\beta$" = estimate_year, 
         "Standard error" = se_year, 
         "Lower bound" = ci_low, 
         "Upper bound" = ci_upper) %>% 
  knitr::kable(digits = 3) %>% 
  kable_styling(full_width = F, bootstrap_options = c("hover", "striped")) %>% 
  footnote(general = "Splines correspond with dates in which 3 major policy changes were inacted. \nEstimates are adjusted for organ donor organization")

Term	β	Standard error	Lower bound	Upper bound
Time (year)	3.088	0.022	3.045	3.131
2012 spline	-0.982	0.041	-1.063	-0.902
2016 spline	0.333	0.150	0.039	0.627
2017 spline	2.465	0.213	2.048	2.883
Note:
Splines correspond with dates in which 3 major policy changes were inacted. Estimates are adjusted for organ donor organization

add_predictions(organ_sp, all_sp_model) %>% 
  ggplot(aes(x = total_days, y = pred, color = opo, type = county)) + 
  geom_line() +
  viridis::scale_color_viridis(discrete = TRUE) + 
  geom_vline(xintercept = 1491, linetype = "dashed") + 
  geom_vline(xintercept = 2799, linetype = "dashed") + 
  geom_vline(xintercept = 3088, linetype = "dashed") + 
  annotate(geom = "text", x = 1720, y = 5, label = "2012") +
  annotate(geom = "text", x = 2600, y = 5, label = "2016") +
  annotate(geom = "text", x = 3300, y = 9, label = "2017") +
  theme(legend.position = "bottom", 
        legend.title = element_blank()) + 
  guides(color = guide_legend(ncol = 2)) + 
  labs(title = "Proportion of organ donors in NY counties", 
       x = "Days since start of data collection",
       y = "% enrolled as donors")

mods <- list(no_sp_model, sp_12_model, sp_16_model, sp_17_model, all_sp_model)

names(mods) <- c("no_sp_model", "sp_12_model", "sp_16_model", "sp_17_model", "all_sp_model")

rmse_fun <- function(model) {
  rmse(model, organ_sp)
}

purrr::map(mods, rmse_fun) %>% 
  as_tibble() %>% 
  gather(key = "model", 
         value = "rmse") %>% 
  mutate(model = ifelse(model == "no_sp_model", "No spline", model), 
         model = ifelse(model == "sp_12_model", "2012 spline", model), 
         model = ifelse(model == "sp_16_model", "2016 spline", model), 
         model = ifelse(model == "sp_17_model", "2017 spline", model), 
         model = ifelse(model == "all_sp_model", "All splines", model)) %>% 
  rename("Model" = model, 
         "RMSE" = rmse) %>% 
  knitr::kable(digits = 3) %>% 
  kable_styling(full_width = F, bootstrap_options = c("hover", "striped"))

Model	RMSE
No spline	1.915
2012 spline	1.902
2016 spline	1.883
2017 spline	1.855
All splines	1.770

options(knitr.kable.NA = '')

anova(no_sp_model, sp_17_model) %>% 
  broom::tidy() %>% 
  select(term, AIC, logLik, statistic, Chi.Df, p.value) %>% 
  mutate(term = ifelse(term == "no_sp_model", "No spline", term), 
         term = ifelse(term == "sp_17_model", "2017 spline", term), 
         p.value = ifelse(p.value < 0.001, "< 0.001", p.value)) %>% 
  knitr::kable(col.names = c("Model", "AIC", "Log-like", "$\\chi^2$", "df", "p-value"), 
               escape = F) %>% 
  kable_styling(full_width = F, bootstrap_options = c("hover", "striped"))

Model	AIC	Log-like	χ²	df	p-value
No spline	30768.94	-15377.47
2017 spline	30306.86	-15145.43	464.084	1	< 0.001

anova(no_sp_model, all_sp_model) %>% 
  broom::tidy() %>% 
  select(term, AIC, logLik, statistic, Chi.Df, p.value) %>% 
  mutate(term = ifelse(term == "no_sp_model", "No spline", term), 
         term = ifelse(term == "all_sp_model", "All splines", term), 
         p.value = ifelse(p.value < 0.001, "< 0.001", p.value)) %>%
  knitr::kable(col.names = c("Model", "AIC", "Log-like", "$\\chi^2$", "df", "p-value"), 
               escape = F) %>% 
  kable_styling(full_width = F, bootstrap_options = c("hover", "striped"))

Model	AIC	Log-like	χ²	df	p-value
No spline	30768.94	-15377.47
All splines	29630.95	-14805.47	1143.996	3	< 0.001

anova(sp_17_model, all_sp_model) %>% 
  broom::tidy() %>% 
  select(term, AIC, logLik, statistic, Chi.Df, p.value) %>% 
  mutate(term = ifelse(term == "sp_17_model", "2017 Spline", term), 
         term = ifelse(term == "all_sp_model", "All splines", term), 
         p.value = ifelse(p.value < 0.001, "< 0.001", p.value)) %>%
  knitr::kable(col.names = c("Model", "AIC", "Log-like", "$\\chi^2$", "df", "p-value"), 
               escape = F) %>% 
  kable_styling(full_width = F, bootstrap_options = c("hover", "striped"))

Model	AIC	Log-like	χ²	df	p-value
2017 Spline	30306.86	-15145.43
All splines	29630.95	-14805.47	679.912	2	< 0.001

# convergence issues with random slope, most likely not needed

# rs <- organ_sp %>% 
#   lmer(eligible_population_enrolled ~ total_days + (1 + total_days | county), data = .) 
# 
# coef(rs)$county %>% 
#   rownames_to_column(var = "county") %>% 
#   mutate(year = total_days * 365)

Mixed Model